Strong coupling theory for driven tunneling and vibrational relaxation 
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We investigate on a unified basis tunneling and vibrational relaxation in driven dissipative multi- 
stable systems described by their N lowest lying unperturbed levels. By use of the discrete variable 
representation we derive a set of coupled non-Markovian master equations. We present analytical 
treatments that describe the dynamics in the regime of strong system-bath coupling. Our findings 
are corroborated by "ab-initio" real-time path integral calculations. 
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Dissipative tunneling in bistable systems finds wide- 
spread applications in many physical and chemical con- 
texts. Usually, the dynamics is described by restricting 
a double-well potential to its lowest tunneling doublet 
of energy eigenstates 0]. This two-level system (2LS) 
describes well the dynamics at low temperature, but be- 
comes increasingly invalid at higher temperatures, when 
higher lying doublets are populated. Then both tunnel- 
ing between the metastable quantum wells and vibra- 
tional relaxation within the wells mix the dynamics. In 
the presence of an external field, e.g., a laser, the differ- 
ent doublets are additionally coupled. Such dissipative 
multilevel systems occur in the context of tunneling of 
magnetization in organic high-spin molecules, i.e., Mni2 
H and Feg ||. Moreover, they mimic tunneling of de- 
fects in metals Q and hydrogen pair transfer in benzoic 
acid [gj. Also, resonant tunneling of the magnetic flux in 
SQUIDS exhibits a multilevel structure Q. An example, 
where a picosecond laser driving accelerates intramolecu- 
lar isomerization in malonaldehyde, is given in |7j . There, 
a pumping mechanism for fast hydrogen transfer ("hy- 
drogen subway") is proposed (see also ||). This complex 
driven dissipative multilevel system has been investigated 
numerically in within a weak coupling approach, and 
in within the real-time path integral tensor multipli- 
cation scheme p"H ] . In absence of driving, the dissipative 
bistable system has been researched by a numerical reac- 
tive flux analytic continuation scheme in p"2| . The only 
analytical approach previous to our work has been at- 
tempted in In these works, however, the dynamics 
restricted to the two lowest tunneling doublets was in- 
vestigated in the weak coupling regime. Because of its 
complexity, even the nondriven multistate dynamics has 
barel y b een studied (l^Jl^]. Common to all prior works 
|§|, [llilJldf| is that the position operator is incorrectly as- 
sumed to be diagonal in the localized (spin) represen- 
tation. This, however, holds true only when the low- 
est tunneling-doublet rules the dynamics, cf. (Q) below. 
Moreover, the harmonic approximation for the two wells 
proposed in Jl3| can be justified only for a large bar- 



rier height where a semiclassical description is already 
applicable Furthermore, the research in |Q is not 
completely based on a full microscopic Hamiltonian. 
The prime objective of this work is the development of a 
first consistent treatment wherein tunneling (TU), driv- 
ing, and vibrational relaxation (VR) are treated analyt- 
ically on a common footing, beyond a Redfield- approach, 
and for an arbitrayry number N of levels. A unified, con- 
sistent treatment becomes possible only if one uses the so 
termed discrete variable representation (DVR) |fiq| . It is 
in this DVR basis that the position operator is diagonal. 
A real-time path integral approach is used to derive a set 
of coupled generalized master equations for the combined 
VR and TU dynamics. The predictions of these equations 
well agree with precise numerical quantum simulations. 
Analytical results for the rate matrix of a JV- level system 
are obtained. 

The model. - To start, we consider the Hamiltonian 
H(t) = Hq + Hext(t) + Hb which accounts for quan- 
tum dissipation and external, generally time-dependent 

2 

control fields. The term H — + V(q) denotes the 
Hamiltonian of the isolated system. We consider a par- 
ticle of mass M moving in the bistable quartic poten- 
tial V(q) = (M 2 U$/64E B )q* - (Afw 2 /4)(j 2 , where E B 
is the barrier height, and loq the frequency of classical 
oscillations around the potential minima at ±qo/2. For 
the isolated system the energy spectrum follows from the 
Schrodinger equation Ho\n) = E n \n), n = 1,2,.... For 
energies well below the barrier, the spectrum consists of 
a ladder of doublets. To illustrate the method, we con- 
sider the case in which only two doublets HAi = E 2 —Ei, 
and ?iA 2 = E 4 — E 3 contribute significantly to the dy- 
namics. They are separated by the energy gap huJo = 
\{Ei + E 3 ) - \{E 2 + Ex) > UAi. We consider then 
the reduction Hq — > H^s, with H^s being the Hamil- 
tonian for the isolated four-level system (4LS). The ex- 
ternal field is characterized by H cxt {t) — —s(t)q, with 
s(t) being a time-dependent field. In the basis of the 
vectors \R X ) = -U|l) + |2)), \Lt) = - |2)), and 



Hh) = ^(|3) + |4)) > \L 2 ) = ^(|3)-|4)), with|^) 
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localized in the left (right) well, the discrete position op- 
erator of the system reads 

i loc = E M^X^-M^X^I) 



A QlCt2 — A a2ai — A l 3 1 f] 2 — A ( g 2i g I = 2v uuj 



(5) 



-b(\L 1 )(R 2 \ + \R 3 ){Li\ 
-\R 1 ){L 2 \-\L 2 ){R 1 \), 



(1) 



with an 



(l|g|2), a 22 = (3|g|4), a 12 = a 21 = 
((l|g|4) + (2|g|3))/2, and & = ((l|g|4) - (2|g|3))/2 « ay. 
Note that, in clear contrast to the 2LS case, it is nondiag- 
onal. For the following analytical treatment, we set b = 0. 
In all results shown in the figures, however, we use 6^0. 
Finally, we model quantum dissipation by an ensemble 
of harmonic oscillators that are bilinearly coupled to the 

2n 



system flf 
with J(u)) 



i.e. 



(7r/2) JV C-/(miUii)6(ui — uJi) being the spec- 



tral density. We assume Ohmic dissipation with an ex- 
ponential cut-off Lu c 3> tuo and a viscous friction strength 
7, such that J {to) — Mju>e~ u '' u ' e . We wish to evaluate 
the probability Pt,(t) := to find at time 

t > to the system in the left well. Here pit) denotes the 
reduced density matrix (RDM) of the system. We bela- 
bor the case in which the particle initially is prepared in 
the lower left state: p(to) = \Li)(Li\, with the bath in 
thermal equilibrium at temperature T = . 

DVR description of vibrational relaxation and 
tunneling. - 

For a harmonic bath, path integral techniques allow one 
to trace out analytically the bath degrees of freedom in 
the eigenbasis of the position operator q (DVR basis). 
Let us start with the introduction of the DVR vectors: 

\ai) = v(\Lx) - u\L 2 )) , |/?i) = vQRx) - u\R 2 )) , 
|oa) = v{u\L x ) + \L 2 )) , \/3 2 ) = uH#i) + \R 2 )) , (2) 
with \ati) (\(3i)) being localized in the left (right) well, 



respectively. Here, v — 1 /vi- 



and u = (an 



A ai )/ai2 = -(022 + A a2 )/ai 2 , and A Qj = -Xp i de- 
note the position eigenvalues: A ai 2 = [—(an + a 22 ) =F 



yj (an — a 22 ) 2 + 4a 2 2 ]/2. Upon the introduction of the 
DVR tunneling matrix elements 

A ail3l = v 2 {A 1 + u 2 A 2 ) , A a2l 3 2 = v 2 (u 2 A 1 + A 2 ) , 

A Qlfe = A Q2/3l ee v 2 u(Ai - A 2 ) , (3) 



being a linear combination of the tunneling splittings Ai 
and A 2 , the 4LS Hamiltonian reads in this DVR basis 

HgV = ~ E 5^ft(k><ftl + l&><«*l) 
+ E + Fpt\0i){0i\) ~ v 2 uhuJ R , (4) 



with F ai = F f3l = u 2 v 2 Hujo, F a2 = Fp 2 = v 2 Tiuj . The 
matrix R = \oti)(a»\ + \a 2 ){ai\ + \Pi)((3 2 \ + |#a)</3i| ac- 
counts for intrawell transitions. It is suggestive to intro- 
duce the DVR intrawell transition elements 



Because A ai p 2 « A ai p l < A a2 p 2 < 2v v&o, different 
time scales determine this VR assisted tunneling dynam- 
ics. For the left well population, one finds Ph(t) = 
Ei^ilpWK)) with the initial RDM p(t ) = = 
u(|ai)(ai| + u 2 \a 2 )(a 2 \ + u\ai)(a 2 \ + u|a 2 )(ai|) being 
nondiagonal. By use of the notation p = ai, a 2 , j3i, (3 2 
for the index of the DVR states, the formally exact path 
integral expression for the diagonal elements of the RDM 
in the DVR basis reads (p^ v := {p\p{t)\v)) 

(*o) , 

with the paths subject to the constraint q(t) = q'(t) = 
X u , and q(t ) = A Mo , q , (t ) = X Uo , with {A„} being the 
position operator eigenvalues. A[q] is the path weight 
in the absence of dissipative forces, while the influence 
functional !F[q, q'] accounts for the bath effects. 

N -level dynamics.- Although derived for a 4LS, these 
equations hold for a finite number N of levels. We switch 
to center of mass r](s) = [q(s) + q'(s)] and relative co- 
ordinates £(s) = [q(s) — g'(s)]. Then, the double path 
integral over the "V-state" paths g(s) or q'{s) is ex- 
pressed as a single path integral over the V 2 states of 
the RDM in the (g, g') plane. The case of the 4LS is 
depicted in Fig. 1. Any such path can be described as 
a sequence of time intervals spent in a diagonal state of 
the RDM ("sojourns") and time intervals between two 
successive visits of the diagonal states ("clusters"). The 
functional T couples different path segments. For a path 
with n transitions at times t\, ..,t n , we introduce the cu- 
mulative off-diagonal charge pj — X)i= i > associated 
to the path derivative £(s) := X)"=i <H S — With 
& = (A Mi - X Ui ) - (A Ms _ 1 - A W4 _J, it yields p n = within 
each cluster of n jumps. Thus, the clusters are "neutral" 
objects, being only weakly interacting at high tempera- 
ture and/or large Ohmic damping, due to the short range 
of the intercluster interaction. This suggests a generaliza- 
tion of the nonintcracting-cluster-approximation (NICA) 
scheme in [T^ ] to our driven situation, with combined 
VR and TU transitions also among non-nearest neigh- 
bors and with an initially nondiagonal RDM. We call this 
the VRTU-NICA. It yields the set of coupled generalized 
master equations (GME) 

puu{t,t Q ) = I u {t,t ) + E / dt'H vli (t,t')p^{t',t ) , (6) 

where all intercluster as well as non-nearest-neighbors 
sojourn-cluster interactions are neglected, while keeping 
fully the intracluster interactions. The VRTU-NICA is 
expected to yield reliable results as long as the friction 
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strength does not exceed level broadening among neigh- 
boring doublets or the temperature is not too low. For 
the case of a single tunneling doublet the DVR basis co- 
incides with the localized basis. Then, the VRTU-NICA 
reduces to the familiar (driven) noninteracting-blip ap- 
proximation for a 2LS [jlj . If more than two levels are in- 
volved, the functions H u ^ and I v are expressed as power 
series in the DVR-Hamilton matrix elements A^. 

Sequential dynamics.- To lowest order one finds 

H v „(t, f) = ^p e -^(*-<') cosfovft - Q'l^t - t')\ , 

L(t,t ) = (5 vot2 - 5 uai )Rep aia2 (t )A aia2 e~ Qa i a 2 {t ~ to) 
x sin[^ QlQ2 (t, t ) - Ql ia2 (t-to)], v^fi. (7) 

The conservation of probability yields for the diagonal 
elements = — X^^u-^M- The driving influence is 
captured by (p Ufi (t,t') = f., dt" [e ^(t") — e u (t")], where 
e v (t) — F v — X u s(t) with F v given below (||). Finally, the 
bath influence is encapsulated in the correlation func- 
tions = Q'py+iQpv with = (A (tl -A y ) 2 Q(i), and 
Q(t) = h Jo°° duj[J{uo)/uj 2 ]{{cosh{(3hu:/2) -cosh[w(^ _ 
it)]} / [sinh((3huj / 2)]) . Hence, to each transition an effec- 
tive friction strength 

aj = (^/qo) 2 a, with a = M7<jg/(27rR) , (8) 

is associated. The prediction for the population Pl(£) 
of the left well is depicted in the inset of Fig. 2, to- 
gether with the results of the quasiadiabatic propaga- 
tor path integral method (QUAPI) adapted to the 
4LS case. For the chosen parameters, higher order co- 
herent paths yield only minor corrections. The inset also 
shows that the dynamics described by ([|) is well approx- 
imated by a Markovian master equation, being indepen- 
dent of the initial off-diagonal preparation, i.e., p vv (t) — 
E^(t)P^{t) > wh ere T ufl (t) = f™ dt' H vft (t,t - t'). 
The explicit time dependence of r,,^ reflects the time 
dependent external forcing. It is appearent that the dy- 
namics is in this regime governed by a single exponen- 
tial decay. To extract this decay rate we observe that 
for high-frequency fields (such as an interdoublet reso- 
nant field), averaging over a driving period is appropri- 
ate. This yields the time-independent rate-matrix 

Kl = ^ j (°°^exphg^(r)]J (^|sin(^)) 

x cosIO^ - F v )t - Q'^(t)] , i/^, (9) 

where C,^ v = A M — A„, with Jq(x) being the zeroth Bessel 
function. The main part of Fig. 2 shows the 4LS av- 
eraged decay rate T av , being the smallest nonzero eigen- 
value of the Markovian rate matrix T^, vs the amplitude 
s of a resonant (f2 = ZJq) driving field s(t) — ssin(Qt). 
Note the perfect agreement between the GME predic- 
tions and those of the QUAPI, together with the charac- 
teristic nonmonotonic behavior. The important issue of 



the contribution of the higher lying states is investigated 
in Fig. 3, where the undriven (s = 0) decay rate is de- 
picted vs the number N of DVR states used to truncate 
the full double-well problem. It is clearly seen that, for 
moderate friction, obeying Aj < 7 < ujq, the truncation 
to the lowest two doublets is adequate even at moder- 
ate temperature fc^T = O.IHujq. This condition implies 
that neighboring doublets do not overlap due to frictional 
broadening. Clearly, the convergence is expected to im- 
prove as the temperature is lowered. We find (not shown) 
that a truncation to a 4LS is adequate for low temper- 
atures and slow driving fields (ft < 0.1 luq). For high- 
frequency fields also higher lying states are involved in 
the dynamics as depicted in the inset of Fig. 3. Hence, 
a reduction of the driven double-well problem to a 4LS 
is problematic in the presence of a strong resonant field 
and moderate temperatures. 

Beyond sequential dynamics. - Upon increasing the tem- 
perature {ksT ~ hu>o), a truncation to a few lev- 
els only starts to be inadequate. Because the effec- 
tive friction strengths aj scale quadratically with £j, 
upon increasing the number N of DVR states involved, 
the iVLS effectively flows to weak coupling. For small 
effective ay, however, the noise action does not sup- 
press long intervals in the off-diagonal states, and the 
higher order paths start to contribute. For 7 < 0.1 loq, 
kgT ~ hujQ and uj c kgT/fi,uj Q , we can approximate 
Q(t) « 2a[ir\t\/K/3 + \n(hf3uj c /2n)} + iua sgni. Now the 
intcrcluster correlations cancel out exactly. The corre- 
sponding averaged Markovian rates read 

x e-i-Ki-V't^Pi/tsfj , (10) 

where pj — X)i=™ &> = 0(1) f° r a vertical (horizontal) 
jump and A tljVj ^ j _ lVj _ 1 is the DVR-Hamilton matrix el- 
ement for the transition from (/ij-i, "j— 1) to [/ij, Vj), e.g. 
for a horizontal jump, K tl]U]4H _ lU] _ 1 := A I/jI/j _ 1 . Since 
( |lC)| ) describes well the high temperature dynamics, this 
rate matrix with s = constitutes the starting point 
for an evaluation of the crossover to the classical regime. 
With s ^ 0, however, a chaotic dynamics may occur. 
Conclusions.- We put forward a real-time path inte- 
gral approach to investigate the interplay between vi- 
brational relaxation and tunneling in driven, dissipative 
multilevel quantum systems. We succeeded in deriving 
a novel GME within the DVR basis which treats tun- 
neling and intrawell dynamics on a common basis: Its 
Markovian approximation in (g) and (|l^) yields novel 
analytical results that well agree with those of the full 
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GME. In contrast to semiclassical imaginary-time rate 
calculations |l5| , we are not limited by the requirement 
of thermal equilibrium at adiabatically varying exter- 
nal fields. Hence, our results provide a powerful tool 
to investigate the crossover from a quantum to a classi- 
cal dynamics beyond the restrictions of the semiclassical 
approximation. Our choice of parameter finely mimics 
the situation of a picosecond laser which accelerates iso- 
merization in malonaldehyde Q| . There a barrier height 
of Eb ~ 1.7 is given, to be contrasted with Eb = 1-4 
here. Also Mni2 and Fes nanomagnets, which have a 
spin ground state of S = 10 (N = 21 levels), provide 

an interesting example where to apply our theory [Q. 
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FIG. 1. The sixteen states of the reduced density matrix of 
a four-level system (4LS) in the discrete variable representa- 
tion. Shown are two typical paths starting and ending in the 
diagonal states O Solid lines indicate vibrational relaxation 
transitions, dotted lines tunneling events. 
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FIG. 2. Averaged tr ansfer r ate of a 4LS vs. the scaled 
(so = hu>o/qo, Qo = WfijMuJo) driving amplitude of a res- 
onant (Q = ZJo = 0.815 u>o) ac-field s sin Qt as obtained from 
the smallest nonzero eigenvalue of the rate matrix (^). The 
three asterisks * denote findings of an exponential fit to the 
QUAPI results. Note the resonant enhancement at finite driv- 
ing strength as the temperature is lowered (To = huio/kB)- 
Here and in the inset we set Eb = lAhuo, 7 = 0.1 uio, 
lo c = 10 loq . Inset: Evolution of the population of the left 
well as predicted by the GME (|) with (Q), by its Marko- 
vian approximation and by the quasi-adiabatic path integral 
method (QUAPI). We choose T = 0.1 T , s = 1.0 s - 
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FIG. 3. Static deacy rate P vs. the number N of levels used 
to truncate the full double-well dynamics. Even at moderate 
temperatures T — 0.1 To the reduction to the two lowest dou- 
blets is appropriate for moderate damping. Here and in the 
inset we set Eb = lAhtoo, uj c = lOtJo- Inset: Averaged trans- 
fer rate T av vs. driving strengths s of a resonant field (fi = ZUo) 
at T = 0.1 Tq. Convergence requires now iV > 4. 
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